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Fluids with competing short-range attractions and long-range repulsions mimic dispersions of charge- 
stabilized colloids that can display equilibrium structures with intermediate range order (IRO), in¬ 
cluding particle clusters. Using simulations and analytical theory, we demonstrate how to detect 
cluster formation in such systems from the static structure factor and elucidate links to macrophase 
separation in purely attractive reference fluids. We find that clusters emerge when the thermal cor¬ 
relation length encoded in the IRO peak of the structure factor exceeds the characteristic lengthscale 
of interparticle repulsions. We also identify qualitative differences between the dynamics of systems 
that form amorphous versus micro-crystalline clusters. 


I. INTRODUCTION 

Complex fluids frequently possess one or more frus¬ 
trating interaction lengthscales that, regardless of ori¬ 
gin, generate micro- to mesoscale structural heterogene¬ 
ity. Archetypical examples include microemulsions [1], 
block copolymers [2, 3], confined fluids [4, 5], and col¬ 
loidal dispersions, including proteins [6-11], wherein the 
surfactant size, block length, pore size, and screened 
electrostatic repulsions set the respective length scales 
of frustration. Despite their contextual differences, all 
exhibit similar transitions between homogeneous fluid 
states and emergent heterogeneous phases with density 
correlations characterized by intermediate range order 
(IRO), typically identified by the presence of a pre-peak 
at low but finite k in the static structure factor S{k) [12]. 

In the case of a pore-glass confined binary fluid sys¬ 
tem [13], the experimental emergence of IRO has been 
rationalized via the behavior of the fluid thermal cor¬ 
relation length ^T, which quantifies the range of corre¬ 
lated concentration fluctuations and the associated IRO 
peak width in 5'(fc). In particular, it was demonstrated 
that the crossover in the temperature-density (T — p) 
plane from dispersed fluid to strong IRO corresponds to 
the conditions at which reaches the pore size, i.e., 
the characteristic frustrating lengthscale. Such condi¬ 
tions enable strong, preferential segregation of the wall- 
attracted species from the other component which, in 
turn, migrates into the pore centers. Additionally, the 
IRO (T — p) crossover conditions corresponded to state 
points close to where the unconfined fluid reference sys¬ 
tem would otherwise exhibit liquid-liquid macrophase 
segregation. 

Here, we extend thermal correlation length concepts 
to a simple model system characterized by IRO: the 
short-range attractive, long-range repulsive (SL) fluid, 
which mimics charge-stabilized colloids with van der 
Waals, depletion, and/or hydrophobic attractions. Var¬ 
ious studies have demonstrated that the long-range re¬ 
pulsive interaction suppresses macrophase separation- 
which would occur for strong short-range attractions 
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alone-in favor of IRO structures including clusters [14- 
18]. However, an ongoing challenge has been to dis¬ 
tinguish between generic IRO (i.e., presence of any 
pre-peak) and clustering specifically, particularly in a 
way accessible to experiments [12, 17]. One such crite¬ 
rion [17] suggests that clustering emerges when the IRO 
peak reaches a magnitude S{k*g^) > 2.7; this bears sim¬ 
ilarity to the empirical Hansen-Verlet single-phase rule 
for tracing macroscopic freezing boundaries in simple 
fluids [19]. 

Section H presents the SL models under consideration 
and the simulation protocol and theoretical methodol¬ 
ogy used to characterize their behaviors. In Section 
HI we propose a new conceptual framework and ac¬ 
curate criterion for clustering: namely, clusters form 
when the thermal correlation length encoded in the 
IRO pre-peak of S{k) exceeds the characteristic length- 
scale of the frustrating interparticle repulsive interac¬ 
tion. We find that this criterion also bolsters previously 
proposed connections between emergent IRO in SL flu¬ 
ids and macroscopic phase separation in corresponding 
reference attractive (RA) models [17] lacking long-range 
repulsions. Finally, we show that the criterion makes 
useful predictions for fluids that form either amorphous 
or micro-crystalline clusters, despite striking qualitative 
differences in the dynamic behaviors of these two types 
of systems. The paper concludes in Section IV with 
a brief summary of our results and their relevance to 
experiment. 


II. METHODS 

Various SL interaction models are known to exhibit 
IRO; here we consider a canonical example given by the 
pairwise potential [14] 


Tsl{x) = 4e(a; 
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xHr 


( 1 ) 


where x = rjd is a, non-dimensionalized particle sep¬ 
aration, d is the measure of particle size, e quantifies 
the attractive strength, and A and respectively char¬ 
acterize the repulsion magnitude and range. We set 
a = 100 in Eq. 1 to mimic archetypical colloids governed 
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by core repulsions with an attraction ranges of 0(1%) 
of the core diameter induced via depletant effects. The 
long-ranged Yukawa tail mimics screened electrostatic 
interactions common to charge-stabilized suspensions. 
The corresponding RA potentials [17] are defined by 
^ra{x) = H{xq — x)ipsl{x)i where H is the Heaviside 
step function and xq is the nearest point for a; > 1 where 
fSL{x) is zero, which eliminates the repulsive tail. 

Model SL fluids defined by Eq. 1 can lose stability 
to micro-crystalline cluster phases at high attraction 
strengths [14], in contrast to many experimental systems 
of interest (e.g., proteins) that do not easily crystallize. 
To study the latter, we also examine a simple ternary 
mixture of SL particles designed to frustrate crystalliza¬ 
tion. The mixture pair potentials are described by 

= 4[e-f(l-2(5,,j)A,](a;-^.“-a;-“)-hA-— 
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where Sij is the Kronecker delta, i,j = —1,0,1 corre¬ 
spond to small, medium (d = 1), and large particles 
respectively, Xij = x — (l/2)(z -|- j)Ad, and pertur¬ 
bative parameter shifts to interaction size and energy, 
Ad = 0.158 and Ag = 0.25, help to thwart crystalliza¬ 
tion and promote mixing, respectively. We use systems 
comprising 20% small, 60% medium, and 20% large par¬ 
ticles. This combination of A^^ and composition repre¬ 
sents a three-component approximation of 10% polydis- 
persity in particle size. 

In examining both models, we set various combi¬ 
nations of the repulsive range and the thermally 
non-dimensionalized repulsive strength pA (where /3 = 
l/ksT and ks is the Boltzmann constant) while varying 
the non-dimensionalized attractive strength /3e. This 
treatment mimics systems for which the short- and 
long-range aspects of constituent interactions are ap¬ 
proximately orthogonal, such as colloids with screening 
lengths set by particle-solvent interactions and attrac¬ 
tions tuned via introduction of depletants [18]. 

To generate equilibrium particle configurations, we 
perform 3D molecular dynamics simulations oi N = 
2960 particles interacting via Eqns. 1 and 2 in the 
NVT ensemble with periodic boundary conditions using 
LAMMPS [20]. Due to the steepness of the repulsion, 
we use an integration time-step of 0.0005, and due to 
the long-range repulsion, we include interactions out to 
a cut-off distance of Xcut = 8.0. For all state points, the 
temperature is fixed at ksT = 1.0 via a Nose-Hoover 
thermostat with time-constant t = 1.0. We calculate 
the structure factor S{k) from simulations by numerical 
Fourier Transform (FT) inversion of the radial distribu¬ 
tion function g{r). To determine whether state points 
are fluid, clustered, or percolating, we calculate cluster 
size distributions (CSDs), which quantify the probabilis¬ 
tic formation of n-particle aggregates, where particles 
are considered part of the same aggregate if their cen¬ 
ters are within the narrow range of the attractive well. 
Similar to other studies [14, 15, 17, 18], a system is con¬ 
sidered clustered with aggregates of preferred size n* 
by the presence of a local maxima in the CSD at n* 
occurring in the range 1 <C n* <C A^, and is consid- 
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FIG. 1. (Color online), (a) Structure factors S{k) for refer¬ 
ence attractive (RA, red dashed) and short-range attractive 
long-range repulsive (SL, blue solid) fluids at packing frac¬ 
tion (p = 0.125 for repulsions with ranges and strengths 
PA. Curves are derived from integral equation theory, where 
the ^R — 10 curves are shown for attraction /?e = 4.35 
and the = 2 curves (shifted vertically) are shown for 
Pe = 4.75. (b,c) S'(fc) curves from (a) replotted to highlight 
fe —^ 0 behaviors, (d) Fourier transforms Poj{k) of the po¬ 
tentials from (a) with ^r = 2 curves shifted vertically, (e) 
RA and SL potentials P^p{r) for the ^r — 2 case. 


ered percolated (at the level of the box) by a CSD peak 
comprised of all particles, i.e., n* ~ N. 

To obtain analytical results for a broader range of po¬ 
tentials, we also derive theoretical thermodynamic and 
pair structure results via the Ornstein-Zernike (OZ) in¬ 
tegral equation relation h(k) = c{k) -|- pc{k)h(k), where 
h{k) = FT[g(r) — 1], c{k) = FT[c(r)], g{r) is the ra¬ 
dial distribution function, c(r) is the direct correlation 
function and p is the number density. The OZ rela¬ 
tion is closed via the Percus-Yevick hard sphere ref¬ 
erence, non-linear optimized random phase approxima¬ 
tion, c(r) « exp[—/3(/3(r)] -1-1- G'(r), where G{r) = 0 
for r > d while for r < c? it is optimized to enforce 
h(r) = — 1 (thus, we approximate Eqn. 1 with a literal 
hard core for r < d) [21]. In carrying out these calcula¬ 
tions, we consider only the Eqn. 1 potential since non¬ 
crystalline states are avoided due to the enforcement of 
homogeneity. This closure yields a spinodal locus at all 
densities, an important feature for the RA cases. 


III. RESULTS AND DISCUSSION 

To begin our discussion, we first consider the behav¬ 
ior of the structure factor S{k) for SL fluids with differ¬ 
ent relative (integrated) repulsive strengths and corre- 
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spending RA systems (see Fig. la-c) as predicted from 
integral equation theory. The two SL fluids exhibit pre¬ 
peaks characteristic of IRO at wavelengths > 0, in¬ 
dicating preferential structuring on microscopic length- 
scales of 2TT/k*g^ « 12.6d and b.Od, respectively. In 
contrast, for the RA fluids lacking long-range repul¬ 
sions, the short-range attractions drive ordering on the 
macroscopic lengthscale, corresponding to the peak at 
k*j^A = 0- Crucially, we see that for the very weak repul¬ 
sive case = 10, pA = 5 x 10“'*) , the S{k) for the SL 
fluid traces the RA curve down to low-A:, supporting the 
conceptual notion of SL fluids as perturbations to un¬ 
derlying RA fluids for which only the principal ordering 
lengthscale has been shifted. 

To understand why one should naturally expect SL 
fluids to aggregate on smaller lengthscales than their 
RA counterparts, we examine in Fig. 1(d) the Fourier 
space analogs of the SL and RA pair potentials, u){k) = 
FT[(^o(?')]) where = H(r — d)ip{r). Viewing the 

potentials in this way makes explicit the idea that struc¬ 
tural oscillations of different lengthscales are weighted 
by the energy profile u}{k), which is evidenced by the 
close reciprocal correspondence between basins in uj{k) 
(Fig. 1(d)) and peaks in S{k) (Fig. 1(a)). This con¬ 
nection can also be made more formal by considering 
microstate configurational energies (see Appendix A). 

Moving beyond the above discussion concerning 
generic IRO, we demonstrate in Fig. 2 that particle 
clustering emerges when the thermal correlation length 
surpasses the characteristic lengthscale of interpar¬ 
ticle repulsion Here, we estimate from the well- 
known S{k) approximation (inverse expansion) near 
k* [21]: 


S{k) = 


S{k*) 

1 + (k — k*Yd?£^’^ 


( 3 ) 


That is a correlation length is evident by consider¬ 
ing the real-space form of Eqn. 3, limr_>.oob(?') — 1] oc 
r“* exp[—r/d^r] cos[rfc*—d], where g{r) is the radial dis¬ 
tribution function, d is a constant, and gives the char¬ 
acteristic decay-length of static correlations while the 
cosine term reflects modulated structure. In practice, 
can be extracted from S{k) by fitting S{k*)/S{k) to 
the form 1 -|- (fc — k*)^d^^^ about k*. 

In Fig. 2(a), we catalog the phase behavior as a func¬ 
tion of attractive strength /3e for various packing frac¬ 
tions 4>. It is evident that for the lower-density iso¬ 
chores, the criterion demarcates when cluster¬ 

ing begins in our polydisperse system, as indicated by 
a characteristic CSD peak with increasing attractions 
(Fig. 2(c)) and reflected by a growing IRO pre-peak in 
S{k) (Fig. 2(d)). As is intuitively expected and seen by 
others [15, 17, 22], for denser isochores like (p = 0.250, 
it is challenging to identify precisely when “clustering” 
begins because the CSD indicates box-wide percolation 
(geometrically merged clusters) even down to relatively 
low /3e. Fig. 2(a) also shows that correlation lengths of 
monodisperse and polydisperse systems coincide upon 
approach to the threshold, where this bound¬ 

ary also approximately identifies where the monodis- 



FIG. 2. (Color online), (a) Symbols show thermal correla¬ 
tion lengths for SL simulations of polydisperse (filled) and 
monodisperse (unfilled) systems with attractive strengths fie 
and packing fractions (j> = 0.050, 0.125, and 0.250. Symbol 
shapes indicate whether the state point is dispersed fluid (tri¬ 
angles), clustered (circle), or percolated (diamond), and the 
horizontal dashed line indicates = Cfi- Solid lines show 
calculated via theory, (b) Phase behavior calculated via the¬ 
ory for potentials from (a), including RA macrophase spin- 
odal (red unfilled squares); SL curves (blue filled squares) 
corresponding to = 2 and = 5; and S{kg]^) = 2.7 
curve (black x). ‘L-l-G’ indicates liquid-gas coexistence, ‘C’ 
indicates clustered phase, and ‘F’ indicates fluid phase, (c) 
Gluster size distributions indicating probability p{n) of n- 
particle cluster formation and (d) S{k) profiles from poly¬ 
disperse simulations at 0 = 0.125. 


perse fluid loses stability with respect to formation of 
micro-crystalline clusters. 

In Fig. 2(b), we also examine phase behaviors for the 
SL and RA fluids derived for a wider Pe — p parame¬ 
ter space via theory, which reveals close correspondence 
between the SL = ^r boundary and the spinodal as¬ 
sociated with RA macrophase separation. Their similar 
shapes (and, in this case, locations), suggest that the 
SL = ^R boundary echoes the RA thermodynamic 
instability, where the frustrating repulsion has erased 
(or highly suppressed) liquid-gas coexistence in favor of 
clustering. (We also include the = 5 curve to demon¬ 
strate the general propagation of the RA spinodal shape 
with increasing /3e.) As a further comparison, the em¬ 
pirical clustering condition S'(/cg^) > 2.7 is also shown. 
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FIG. 3. (Color online), (a) Phase diagrams calculated via 
theory, comprising RA macrophase spinodals (unfilled red 
symbols) and SL curves (filled blue symbols) for 

= 10 and two repulsive strengths PA. (b) RA spinodals 
and curves along which S{kg]^) = 2.7 (filled black symbols) 
for same systems as in (a), (c) Phase diagram calculated 
via theory comprising RA macrophase spinodal (unfilled red 
triangles); SL curves corresponding to macrophase spinodal 
at low (j) (right-pointing blue triangles) and = 2 at 

high (j) (left-pointing blue triangles); and disorder line (pur¬ 
ple squares) in the fluid region (see text). ‘L-|-G’ indicates 
liquid-gas coexistence, ‘C’ indicates clustered phase, and ‘F’ 
indicates fluid phase. 


While it lies within similar proximity to the RA spin¬ 
odal, it possesses a noticeably different, shallower con¬ 
tour. 

To elucidate deeper connections between the contours 
in Fig. 2(b), we explore in Fig. 3 whether the 
and RA spinodal boundaries truly converge for ultra- 
weak repulsions, which might be expected if the latter 
can be considered a natural weak-repulsion limit of the 
former. In Fig. 3(a-b), we examine two potentials with 
different repulsive strengths: for f3A = 1 x 10“^, the re¬ 
pulsion is evidently “strong” and there is no overlap be¬ 
tween the = ^R and RA spinodal boundaries (note: 
this highlights that these boundaries do not generally 
overlap as in Fig. 2(b)). However, as repulsion strength 
is lowered to pA < 1 x 10“®, the two curves collapse 
and become truly indistinguishable, reflecting a deep 
SL-RA connection. In Fig. 3(b), we also show corre¬ 
sponding S{kg^) = 2.7 curves. Clear discrepancies in 
shape are apparent when comparing the RA spinodals 
and the <S'(fcg^) = 2.7 boundaries, and the two types of 
curves increasingly move apart as PA is reduced. 

To further generalize the connection of the RA spin¬ 
odal to the phase behaviors of SL systems, we con¬ 
sider in Fig. 3(c) a less long-ranged weak repulsion 



(a) Polydisperse, (3e=5.2 



(b) Monodisperse, p£=5.2 


FIG. 4. (Color online). Cluster phase simulation snap¬ 
shots of polydisperse (a) and monodisperse (b) systems at 
p = 0.125 with attractive strength /3e = 5.2 and repulsions 
defined by = 2 and PA — 0.20. Particles comprising a 
single cluster (determined at time t) are rendered opaque in 
their positions at times t (left) and t' — t + At (right). The 
lag time is At = 25Td, where Td = d? /D is the characteristic 
time for d = 1 particles to diffuse and D is the long-time 
bulk diffusion coefficient determined via mean-squared dis¬ 
placements. Colors correspond to small, medium (d = 1), 
and large particles, which are shaded yellow, red, and blue, 
respectively. Visualizations created with VMD [23]. 


{pR = 2,PA = 5 X 10“^), which exhibits intriguing 
properties: a true SL spinodal separation occurs for 
(j) < 0.09, while for higher volume fractions there is a 
= Cr clustering boundary. The low-density fluid 
also exhibits a disorder line, below which the IRO peak 
is present and above which the IRO peak transitions to 
a = 0 peak. The intimate correspondence between 
the SL boundaries and the RA spinodal further reflects 
that the condition = pR reflects a muted thermo¬ 
dynamic instability, which for very weak repulsions can 
also emerge within the SL fluid itself. 

Finally, we consider the morphologies and lifetimes of 
the clusters that form in polydisperse and monodisperse 
SL systems. Clusters in the former exhibit amorphous 
and irregular shapes, as exemplified by the simulation 
snapshots in Fig. 4(a), which correspond to the sys¬ 
tem in Fig. 2 at conditions slightly above the clustering 
transition. Here, it is evident based on the time-lag 
snapshots that the clusters are transient and contin¬ 
uously redistribute particles to create new clusters at 
the expense of others. By significantly increasing the 
attractive strength /3e, one can eventually observe ar¬ 
rested, percolating, amorphous gels as exemplified by 
the simulation snapshots for p = 0.125 systems in Fig. 
5. Interestingly, our model gels may be thermoreversible 
with no local crystallinity, possibly providing a simpler 
alternative to valence-limited gel-formers [24]. Ther- 
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moreversibility is highly desired to facilitate fabrication 
of massively reconfigurable, reversible materials. 

In contrast, monodisperse systems at similar attrac¬ 
tion strengths can undergo highly regular clustering 
via local crystallization, as exemplified in Figure 4(b). 
While the crystalline nature of such simulated clusters 
has been observed previously by others [14, 15, 18], we 
do note that the relatively weaker repulsion examined 
here drives formation of much larger clusters that are 
more obviously crystalline in nature. The crystalline 
clusters are relatively static objects once formed, as 
demonstrated by the time-lag snapshots, in direct con¬ 
trast to the amorphous clusters. 



(a) «i=0.125, pf;=5.50 



(b) «i=0.125, pe=5.80 


(c) ^=0.125, pe=6.30 (d) ^=0.125, pe=6.80 

FIG. 5. (Color online). Cluster phase simulation snap¬ 
shots of polydisperse systems at = 0.125 with various at¬ 
tractive strengths fit and repulsions defined by = 2 and 
pA = 0.20. In all snapshots, particles comprising a sin¬ 
gle cluster (determined at time t) are rendered opaque in 
their positions at time t. For cases (a) and (b) that are not 
gelled, the same particles are also shown in their positions 
aXt'=t-\- At. The lag time in (a) and (b) At = 25rd, where 
Td = dP/D is the characteristic time for d = 1 particles 
to diffuse and D is the long-time bulk diffusion coefficient 
determined via mean-squared displacements. For cases (c) 
and (d), the configurations are dynamically arrested and 
cannot be practically measured within the timescale of sim¬ 
ulations. Colors correspond to small, medium (d = 1), and 
large particles, which are shaded yellow, red, and blue, re¬ 
spectively. Visualizations created with VMD [23]. 


IV. CONCLUDING REMARKS 

In closing, we have presented a new framework for 
understanding and detecting cluster phases in SL fluids 
based on the thermal correlation length px- This frame¬ 
work should prove useful for probing micro-structural 
transitions in diverse systems governed by frustrated 
interactions, e.g., lattice spin models with opposing 
nearest-neighbor and higher-order couplings. We have 
also presented the first non-microcrystallizing SL fluid, 
which exhibits amorphous transient clusters; this should 
prove useful for examining the (zeroth order) physics of 
real dispersions known to be resistant to crystallization, 
e.g., proteins. 

Finally, we remark that the Pt = clustering cri¬ 
terion can be implemented in experiments provided 
that, in addition to extracting px from an S{k) pro¬ 
file (described earlier), one can also obtain a reason¬ 
able measure of the repulsive lengthscale between par¬ 
ticles pR. For systems accurately described by simple 
screening models, Pr can be directly estimated. Oth¬ 
erwise, one can first obtain the r-space total correla¬ 
tion function h{r) via an inverse FT of S{k). Like¬ 
wise, one can calculate the direct correlation func¬ 
tion c{k) = p~^ — [p5'(A:)]“^ and then obtain its r- 
space equivalent c(r) = FT~^[c(fc)], which provides in¬ 
formation about the interparticle interactions because 
limr_).oo c(r) « p{r) [21]. By plotting ln{jr/i(r)j} and 
ln{jrc(r)j} versus r (where jccj is the absolute value of 
x) and comparing their (negative) slopes, one directly 
compares the range of interparticle correlations (as cap¬ 
tured by ^t) and the characteristic range of the interpar¬ 
ticle interactions, respectively. Thus, given an S{k) pro¬ 
file exhibiting an IRO peak, if ln{jr/i(r)j} decays more 
slowly than ln{jrc(r)j}, then the ^x associated with IRO 
exceeds the characteristic (repulsive) lengthscale 
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APPENDIX A: FOURIER-SPACE 
CONNECTIONS BETWEEN io{k) AND S{k) 

An N particle configuration [r^] that does not violate 
the hard core constraint is weighted according to the 
Boltzmann factor exp[—/3U([ri])], where: 

1 ^ 

^([u]) = x ‘^odu-i-j!) (Ai) 

is the total potential energy due to the non-hard-core 
portion of the pair potential poir). Eqn. Al can be 
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recast using the definition of the 3D dirac delta function ^(x): 


D([r,]) = 1'^ / - Ri)¥^o(|Ri - R2|)'5(r, - R 2 ) (A2) 


Since Eqn. A2 is a convolution with respect to Ri 
and R 2 , it can be recast as a single integral in Fourier 
space using the Fourier transformed potential a;(fc) = 
FT[(^o(r)]: 

1 ^ 1 /• 

0([r,]) ^ (A3) 

2 ( 2 ^)' J 

Moving the sum inside the integral in Fqn. A3 and 
using the definition of the non-ensemble averaged total 
correlation function, 

N 

h{k-, [r,]) = {pNy^ exp[-ik • (r, - r^)] (A4) 


one can subsequently write 


^^([ d ]) = J dku;ik)h{k-, [r,]) (A5) 

which makes explicit the role uj{k) plays in favoring [r^] 
states possessing certain oscillatory structural correla¬ 
tions. Namely, any thermodynamically favorable config¬ 
uration [r*], as weighted by exp[—/3D([ri])], is captured 
by the equilibrium average total correlation function 
h{k) « h{k] [r*]). In turn, ijj{k) sets the energetic “pref¬ 
erence” for configurations structured at certain wave¬ 
lengths fc, which appear as peaks in the structure factor 
since S{k) = 1 + ph{k). 
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